# ======= analyse data Accuracy II, experiment 1 ======= #
# Folco Panizza, 2022/10/03

# load libraries ##############
library(data.table)           # extention of data.frame
library(here)                 # find your project's files, based on the current working directory at the time when the package is loaded
setDTthreads(threads = 0)     # number of threads to be used for data.table parallelised operations
library(rms)                  # a collection of functions for binary and ordinal logistic regression models
library(multcomp)             # simultaneous tests and confidence intervals for general linear hypotheses in parametric models https://stats.idre.ucla.edu/r/faq/how-can-i-test-contrasts-in-r/
library(ordinal)              # Implementation of cumulative link (mixed) models also known as ordered regression models
library(emmeans)              # Obtain estimated marginal means (EMMs) for many linear, generalized linear, and mixed models
library(lmerTest)             # Tests in linear mixed effects models
library(boot)                 # functions and datasets for bootstrapping


# graphical libraries #########
graphics <- c('ggplot2',      # a system for 'declaratively' creating graphics, based on "The Grammar of Graphics".
              'ggforce',      # ggplot2 extension that provides facilities for composing specialised plots 
              'ggrepel',      # provides text and label geoms for 'ggplot2' that help to avoid overlapping text labels.
              # 'GGally',       # functions for pairwise plot matrices
              'scales',       # graphical scales map data to aesthetics, and provide methods for automatically determining breaks and labels for axes and legends.
              'ggpubr',       # provides some easy-to-use functions for creating and customizing 'ggplot2'-based publication ready plots
              'extrafont',    # tools to using fonts other than the standard PostScript fonts. 
              'ggthemes',     # extra themes for ggplot2
              'RColorBrewer') # Provides color schemes for maps (and other graphics)

invisible(lapply(graphics, library, character.only = TRUE)); rm(graphics)

# set default theme for plots
theme_set(theme_classic() + theme(text=element_text(size=24), panel.grid.major = element_line(colour="grey90"))) #family="Trebuchet MS", 


# custom functions ============

backup_options <- options(); options(digits=3) # print 3 meaningful digits
# options(backup_options)

# summary with 95% confidence interval
summary2 <- function(x) {m<-lm(x ~ 1); cbind(coef(m),confint(m))} # c(Estimate=mean(x), `Std. Error`=sd(x)/sqrt(length(x)))

# summary for multiple contrasts
comparison.summary <- function(comp, fwe="none"){
  table.test <- as.data.frame(summary(comp, test = adjusted(fwe))$test[c("coefficients","tstat","pvalues")])
  table.ci   <- confint(comp, calpha = if(fwe!="none"){adjusted_calpha()}else{univariate_calpha()})$confint[,-1]
  if(is.null(dim(table.ci))) table.ci <- t(table.ci)
  table.fit  <- cbind(table.test,table.ci)[,c("coefficients","lwr","upr","tstat","pvalues")] 
  setDT(table.fit, keep.rownames = T); setnames(table.fit,"pvalues","p.value")
  table.fit[,`   `:=fcase(p.value>.1,"   ",p.value>.05,".  ",p.value>.01,"*  ",p.value>.001,"** ",default="***")]
  table.fit[,colnames(table.fit)[2:6] := lapply(.SD, round, 3), .SDcols = 2:6]
  return(table.fit[])}

# load datasets ===============
load(here("Data","experiment 1","questionnaire.Rdata")) # data from qualtrics
load(here("Data","experiment 1","search-data.Rdata"))   # data from lab.js

# add post timing to Qualtrics data
post_timing = JavaScript[, .(post_evaluation = unique(active_time)/1000), .(ResponseId)]
Qualtrics = merge(Qualtrics, post_timing, by = "ResponseId"); rm(post_timing)

# Exclusion criteria ==========
# We will exclude from the analyses only participants.
mobile = JavaScript[,unique(mobile_user),ResponseId][(V1),ResponseId] # who are using a mobile device
# length(mobile) # 4
Qualtrics  = Qualtrics [!(ResponseId%in%mobile)]
JavaScript = JavaScript[!(ResponseId%in%mobile)]

# save(Qualtrics, file = "C:/Users/folco/Google Drive/Peritia 2020 - Milan Team/Research Projects/Accuracy II/Meta/study2exp1.Rdata")

# rsdnce = Qualtrics[residence!="United Kingdom", residence] # who are not U.K. residents
# length(rsdnce) # 0

# Qualtrics[Clarity=="Yes", cat(paste(1:.N,Clarity_5_TEXT, collapse = "\n"))]     # Comments
# Qualtrics[Technical=="Yes", cat(paste(1:.N,Technical_5_TEXT, collapse = "\n"))] # Technical issues few people reported technical issues outside a clear lag in Qualtrics, non were "invalidating"
# external search: comments
# Qualtrics[Where %like% "Other" , cat(paste("-",unique(Where_6_TEXT), collapse = "\n"))]
# Qualtrics[WhyNot %like% "Other", cat(paste("-",unique(WhyNot_5_TEXT), collapse = "\n"))]

# Descriptive analyses ========

Qualtrics[,{d <- table(incentive); print(d); chisq.test(d)}] # design balance
# Qualtrics[,{d <- table(tag); print(d); chisq.test(d)}] # by-post, p>.999
# control monetary     self   social 
#     992     1003     1004     1000 
# Chi-squared test for given probabilities
# X-squared = 0.09, df = 3, p-value = 1

# timing
Qualtrics[, .(medianDuration = median(time_taken/60)), incentive] # of the whole study (in minutes) # `Duration (in seconds)` variable: same but from Qualtrics
#    incentive medianDuration
# 1:      self           4.09
# 2:  monetary           4.39
# 3:    social           4.16
# 4:   control           4.00

Qualtrics[, range(time_taken/60)]
# 0.514 98.792

Qualtrics[, as.list(summary(post_evaluation)), incentive]
#    incentive Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1:      self 3.68    22.4   34.6 58.5    64.5 1173
# 2:  monetary 3.62    25.4   41.8 67.2    81.6  695
# 3:    social 5.68    21.3   34.5 59.5    64.1  847
# 4:   control 3.49    21.4   32.2 53.5    56.0 1471

# pounds/hour
Qualtrics[, median((.70 + fifelse(incentive=="monetary" & correct==1, .70, 0))/time_taken*3600)]
# 11.2

# accuracy
Qualtrics[, {s <- summary2(correct)       ; split(s, colnames(s))}, incentive]
#    incentive       2.5 % 97.5 %
# 1:      self 0.651 0.622  0.681
# 2:  monetary 0.705 0.677  0.733
# 3:    social 0.632 0.602  0.662
# 4:   control 0.612 0.582  0.642

Qualtrics[, {s <- summary2(accuracy_score); split(s, colnames(s))}, incentive]
#    incentive      2.5 % 97.5 %
# 1:      self 4.01  3.92   4.10
# 2:  monetary 4.22  4.13   4.31
# 3:    social 3.95  3.87   4.04
# 4:   control 3.87  3.78   3.95

# demographics
Qualtrics[, table(droplevels(Sex))/.N] # gender
#  Female              Male Prefer not to say 
# 0.62666           0.37109           0.00225 
Qualtrics[, .(M=mean(age, na.rm = T), SD=sd(age, na.rm = T), `NA`=sum(is.na(age)))] # age
#    M   SD NA
# 37.5 13.4  6
Qualtrics[education%like%"grad|doc", .N/nrow(Qualtrics)*100] # bachelor education or higher
# 58.7

# Preregistered analyses ======

contrasts.lrm <- c("`incentive=monetary` = 0", # 3. Main effect of monetary incentives on accuracy.
                   "`incentive=self` = 0",     # 5. Main effect of accuracy comparison on accuracy.
                   "`incentive=social` = 0",   # 4. Main effect of social responsibility on accuracy.
                   # comparison between conditions
                   "`incentive=monetary` - `incentive=self` = 0",
                   "`incentive=monetary` - `incentive=social` = 0",
                   "`incentive=self` - `incentive=social` = 0")

contrasts.glm <- c("incentivemonetary = 0","incentiveself = 0","incentivesocial = 0","incentivemonetary - incentiveself = 0","incentivemonetary - incentivesocial = 0","incentiveself - incentivesocial = 0")

# accuracy
fit  <- lrm(accuracy_score~incentive, Qualtrics)

comp <- glht(fit, linfct = contrasts.lrm) #[1:3] research question 1 #[4:5] research question 3
comparison.summary(comp, "fdr")
#                                           rn coefficients    lwr   upr tstat p.value    
# 1:                        incentive=monetary        0.495  0.290 0.699  6.21   0.000 ***
# 2:                            incentive=self        0.201 -0.002 0.404  2.54   0.016 *  
# 3:                          incentive=social        0.114 -0.087 0.315  1.45   0.175    
# 4:   `incentive=monetary` - `incentive=self`        0.294  0.088 0.500  3.66   0.000 ***
# 5: `incentive=monetary` - `incentive=social`        0.381  0.176 0.585  4.79   0.000 ***
# 6:     `incentive=self` - `incentive=social`        0.087 -0.116 0.290  1.10   0.270    

# +4.7% +1.1% +7.8% # estimated increase in proportion of correct answers, self condition
logodds = data.frame(Estimate = coef(fit), Std.Error = sqrt(diag(vcov(fit)))); logodds = with(logodds[c("y>=4","incentive=self"),],Estimate + cbind(Std.Error)%*%qnorm(c(.5, .025, .975)))
diag(plogis(outer(logodds[1,], logodds[2,], `+`)) - matrix(rep(plogis(logodds[1,]),3), 3, byrow = T))

# accuracy plot ++++++++++++++
sina_plot = Qualtrics[, boot(.SD, function(x, i) {x[i, c(accuracy = mean(accuracy))]}, R = 1000)$t, by = incentive]
sina_plot[, incentive := factor(incentive, levels = c("control","social","self","monetary"), labels = c("control","inform\nothers","peer\nperformance","monetary\nbonus"))]

stat_label  = data.frame(incentive=c("control","control","inform\nothers","peer\nperformance"), V1 = c(.665,.705,.715,.725), label=c("*","***","***","***"))

# cairo_pdf(here("Figures","experiment 1","incentives-score.pdf"), width = 12, height = 5)
png(here("Figures","experiment 1","incentives-score.png"), width = 12, height = 4.5, units = "in", res = 216)
ggplot(sina_plot, aes(x = incentive, y = V1, fill = incentive, colour = incentive)) + 
  theme(legend.position = "top", legend.key.width = unit(4, "lines"), legend.text=element_text(margin=margin(l=-25)),
        axis.line.y = element_blank(), axis.text.x = element_text(size = 20), axis.ticks.x = element_line(size = .75),
        panel.grid.major.x = element_line(size=.1, color="grey65"), panel.grid.minor.x = element_line(size=.1, color="grey65"), aspect.ratio = 2/7) + 
  geom_sina(alpha=0.05, scale=F, method="density", maxwidth = .9, size = 4, shape = 19) + 
  geom_hline(yintercept = 0.5, linetype = 2, size = 1, color="grey25") +
  stat_compare_means(data = sina_plot, label.y = c(.65,.69,.70,.71), tip.length = 0, comparisons = list(c("peer\nperformance", "control"),c("monetary\nbonus", "control"),c("monetary\nbonus", "inform\nothers"),c("monetary\nbonus", "peer\nperformance")), bracket.size = 1, size = 7, symnum.args = list(cutpoints = c(0.000,.5,1.000), symbols = c("",""))) +
  geom_text(data = stat_label, label=stat_label$label, color = "black", size = 12, vjust = 1, hjust = 0) +
  annotate("text", x = 2.5, y = .5, label = "random response", angle = 90, vjust = -1, color = "grey25", size = 5.5) +
  coord_flip() +
  scale_fill_manual(values=c("#D55E00", "#00D199", "#009E73","#006B4F"), guide = FALSE, name = element_blank()) +
  scale_colour_manual(values=c("#D55E00", "#00D199", "#009E73","#006B4F"), guide = FALSE, name = element_blank()) +
  scale_x_discrete(breaks = NULL, name = NULL, limits = rev(levels(sina_plot$incentive))) + 
  scale_y_continuous(name = "accuracy score", limits = c(.4, 1), breaks = c(.4, .6, .8, 1), labels = function(x) x*5+1, minor_breaks = c(.7, .9), expand = expansion(mult = c(.04, .04))) +
  geom_text(data = data.frame(incentive = levels(sina_plot$incentive), V1 = c(.53, .53, .55, .59)), label = levels(sina_plot$incentive), size = 4.5, fontface = 'bold') +
  geom_label_repel(data = sina_plot[,mean(V1), incentive], mapping = aes(label = round(V1*5+1, 2)), fontface = 'bold', color = 'black', fill = 'white', max.iter = 3e2, box.padding = 1, point.padding = 0, segment.color = NA, force = 0)
dev.off()

rm(sina_plot, stat_label)


# performance by technique
Qualtrics[lateral==F, {m <- summary2(accuracy_score); split(m,colnames(m))}, incentive][order(incentive)]
Qualtrics[restraint==T, {m <- summary2(accuracy_score); split(m,colnames(m))}, incentive][order(incentive)]

# technique use by condition
# Qualtrics[,table(incentive,fcase(lateral==F,"non-adopter",restraint==T,"adopters",default = "lateral-only"))]
# Qualtrics[,prop.table(table(incentive,fcase(lateral==F,"non-adopter",restraint==T,"adopters",default = "lateral-only")),1)]
# incentive  adopters lateral-only non-adopter
#  control     0.0524       0.0353      0.9123
#  monetary    0.0927       0.0588      0.8485
#  self        0.0548       0.0378      0.9074
#  social      0.0670       0.0320      0.9010

fit <- lrm(accuracy_score ~ lateral + incentive, Qualtrics)
comp <- glht(fit, linfct = c("`lateral=TRUE`=0")) #[1:3] research question 1 #[4:5] research question 3
comparison.summary(comp, "fdr")
#              rn coefficients   lwr   upr tstat p.value    
# 1: lateral=TRUE        0.431 0.252 0.609  4.73       0 ***

fit <- lrm(accuracy_score ~ restraint + incentive, Qualtrics)
comp <- glht(fit, linfct = c("`restraint=TRUE`=0")) #[1:3] research question 1 #[4:5] research question 3
comparison.summary(comp, "fdr")
#                rn coefficients  lwr   upr tstat p.value    
# 1: restraint=TRUE        0.542 0.32 0.764  4.78       0 ***

logodds = data.frame(Estimate = coef(fit), Std.Error = sqrt(diag(vcov(fit))))
# logodds = with(logodds,Estimate + cbind(Std.Error)%*%qnorm(c(.5, .025, .975)))
# 
# plogis(logodds[3,]) # model predictions for the control condition
# # 60.5% [57.8%, 63.2%] # estimate
# # 61.2% [58.1%, 64.2%] # estimate glm
# # 61.2% [58.2%, 64.2%] # lm
# 
# diag(plogis(outer(logodds[3,], logodds[8,], `+`))) # model predictions for the inform others condition
# # 63.2% [56.8%, 69.2%] # estimate
# # 63.2% [55.8%, 70.0%] # estimate glm
# # 63.2% [60.2%, 66.2%] # +/- 1.96*SE
# 
# diag(plogis(outer(logodds[3,], logodds[7,], `+`))) # model predictions for the challenge skills condition
# # 65.2% [58.9%, 71.0%] # estimate
# # 65.1% [57.8%, 71.8%] # estimate glm
# # 65.1% [62.2%, 68.1%] # +/- 1.96*SE
# 
# diag(plogis(outer(logodds[3,], logodds[6,], `+`))) # model predictions for the monetary bonus condition
# # 71.6% [65.8%, 76.7%] # estimate
# # 70.5% [63.6%, 76.6%] # estimate glm
# # 70.5% [67.7%, 73.3%] # +/- 1.96*SE

logodds = with(logodds[c("y>=4","restraint=TRUE"),],Estimate + cbind(Std.Error)%*%qnorm(c(.5, .025, .975)))
# logodds = with(logodds[c("y>=6","restraint=TRUE"),],Estimate + cbind(Std.Error)%*%qnorm(c(.5, .025, .975)))
diag(plogis(outer(logodds[1,], logodds[2,], `+`)) - matrix(rep(plogis(logodds[1,]),3), 3, byrow = T)) # diag()
# proportion of 4,5 or 6 answers:
# +12.1% [+7.6%,+15.6%]

# estimate using logistic regression (different statistical assumptions)
# fit <- glm(correct ~ restraint + incentive, Qualtrics, family = "binomial")
# logodds = as.data.frame(coef(summary(fit))[,1:2])
# logodds = with(logodds[c("(Intercept)","restraintTRUE"),],Estimate + cbind(`Std. Error`)%*%qnorm(c(.5, .025, .975)))
# diag(plogis(outer(logodds[1,], logodds[2,], `+`)) - matrix(rep(plogis(logodds[1,]),3), 3, byrow = T))
# proportion of 4,5 or 6 answers:
# +10.2% [+4.1%,+14.9%]



fit <- glm(lateral ~ incentive, Qualtrics, family = "binomial")
#                                     rn coefficients    lwr   upr  tstat p.value    
# 1:                   incentivemonetary        0.619  0.253 0.986  4.342   0.000 ***
# 2:                       incentiveself        0.060 -0.341 0.461  0.384   0.701    
# 3:                     incentivesocial        0.134 -0.262 0.530  0.866   0.580    
# 4:   incentivemonetary - incentiveself        0.559  0.200 0.919  3.995   0.000 ***
# 5: incentivemonetary - incentivesocial        0.486  0.132 0.839  3.528   0.001 ***
# 6:     incentiveself - incentivesocial       -0.074 -0.463 0.316 -0.484   0.701    

fit <- glm(restraint ~ incentive, Qualtrics, family = "binomial")
#                                     rn coefficients    lwr   upr  tstat p.value    
# 1:                   incentivemonetary        0.614  0.154 1.074  3.423   0.004 ** 
# 2:                       incentiveself        0.047 -0.464 0.557  0.234   0.815    
# 3:                     incentivesocial        0.261 -0.228 0.750  1.370   0.256    
# 4:   incentivemonetary - incentiveself        0.567  0.115 1.020  3.217   0.004 ** 
# 5: incentivemonetary - incentivesocial        0.353 -0.075 0.781  2.115   0.069 .  
# 6:     incentiveself - incentivesocial       -0.214 -0.696 0.267 -1.142   0.304    

comp <- glht(fit, linfct = contrasts.glm)
comparison.summary(comp, "fdr")



# Exploratory analyses ========

# correct guessing
fit <- glm(correct ~ incentive, Qualtrics, family = "binomial")

# estimated increase in the proportion of correct answers, self condition
logodds = as.data.frame(coef(summary(fit))[,1:2]); logodds = with(logodds[c("(Intercept)","incentiveself"),],Estimate + cbind(`Std. Error`)%*%qnorm(c(.5, .025, .975)))
diag(plogis(outer(logodds[1,], logodds[2,], `+`)) - matrix(rep(plogis(logodds[1,]),3), 3, byrow = T))
# +4.0% -0.3% +7.6%
# +4.7% +1.1% +7.8% # estimate from `lrm` function

comp <- glht(fit, linfct = contrasts.glm)
comparison.summary(comp, "fdr")
#                                     rn coefficients    lwr   upr tstat p.value    
# 1:                   incentivemonetary        0.415  0.171 0.660 4.369   0.000 ***
# 2:                       incentiveself        0.170 -0.069 0.409 1.829   0.101    
# 3:                     incentivesocial        0.086 -0.152 0.323 0.925   0.365    
# 4:   incentivemonetary - incentiveself        0.245 -0.001 0.492 2.562   0.021 *  
# 5: incentivemonetary - incentivesocial        0.330  0.085 0.575 3.459   0.002 ** 
# 6:     incentiveself - incentivesocial        0.084 -0.155 0.324 0.905   0.365    
comparison.summary(comp, "none")
#                                     rn coefficients    lwr   upr tstat p.value    
# 1:                   incentivemonetary        0.415  0.229 0.602 4.369   0.000 ***
# 2:                       incentiveself        0.170 -0.012 0.352 1.829   0.067 .  
# 3:                     incentivesocial        0.086 -0.096 0.267 0.925   0.355    
# 4:   incentivemonetary - incentiveself        0.245  0.058 0.433 2.562   0.010 *  
# 5: incentivemonetary - incentivesocial        0.330  0.143 0.517 3.459   0.001 ***
# 6:     incentiveself - incentivesocial        0.084 -0.098 0.267 0.905   0.365    


# response times
fit <- lm(log(post_evaluation)~incentive, Qualtrics)
# fit <- lm(frank(post_evaluation)~incentive, Qualtrics)
# fit <- lm(log(log(post_evaluation))~incentive, Qualtrics)
comp <- glht(fit, linfct = contrasts.glm)
comparison.summary(comp, "fdr")
#                                     rn coefficients    lwr   upr tstat p.value    
# 1:                   incentivemonetary        0.249  0.155 0.342 6.821   0.000 ***
# 2:                       incentiveself        0.083 -0.010 0.177 2.289   0.033 *  
# 3:                     incentivesocial        0.068 -0.026 0.162 1.867   0.074 .  
# 4:   incentivemonetary - incentiveself        0.165  0.072 0.259 4.547   0.000 ***
# 5: incentivemonetary - incentivesocial        0.181  0.087 0.274 4.963   0.000 ***
# 6:     incentiveself - incentivesocial        0.015 -0.078 0.109 0.421   0.674    



# To test for robustness, we will test if main results hold excluding participants who:
# - failed the attention check measures;
Qualtrics[attentionchecks == 2, .N] # -2: failed attention check; 2: passed check 

fit  <- lrm(accuracy_score~incentive, Qualtrics, subset = attentionchecks == 2)
comp <- glht(fit, linfct = contrasts.lrm)
#                                           rn coefficients    lwr   upr tstat p.value    
# 1:                        incentive=monetary        0.502  0.297 0.708  6.28   0.000 ***
# 2:                            incentive=self        0.203  0.000 0.407  2.57   0.015 *  
# 3:                          incentive=social        0.111 -0.091 0.313  1.41   0.189    
# 4:   `incentive=monetary` - `incentive=self`        0.299  0.092 0.506  3.71   0.000 ***
# 5: `incentive=monetary` - `incentive=social`        0.391  0.186 0.596  4.90   0.000 ***
# 6:     `incentive=self` - `incentive=social`        0.092 -0.111 0.296  1.17   0.244 

fit <- glm(lateral ~ incentive, Qualtrics, family = "binomial", subset = attentionchecks == 2)
comp <- glht(fit, linfct = contrasts.glm)
#                                     rn coefficients    lwr   upr  tstat p.value    
# 1:                   incentivemonetary        0.618  0.252 0.984  4.328   0.000 ***
# 2:                       incentiveself        0.045 -0.358 0.447  0.286   0.775    
# 3:                     incentivesocial        0.131 -0.265 0.528  0.851   0.592    
# 4:   incentivemonetary - incentiveself        0.573  0.213 0.934  4.079   0.000 ***
# 5: incentivemonetary - incentivesocial        0.486  0.133 0.840  3.530   0.001 ***
# 6:     incentiveself - incentivesocial       -0.087 -0.477 0.304 -0.569   0.683    

fit <- glm(restraint ~ incentive, Qualtrics, family = "binomial", subset = attentionchecks == 2)
comp <- glht(fit, linfct = contrasts.glm)
#                                     rn coefficients    lwr   upr  tstat p.value    
# 1:                   incentivemonetary        0.612  0.151 1.073  3.412   0.003 ** 
# 2:                       incentiveself        0.024 -0.489 0.537  0.119   0.905    
# 3:                     incentivesocial        0.259 -0.231 0.748  1.358   0.255    
# 4:   incentivemonetary - incentiveself        0.588  0.132 1.044  3.317   0.003 ** 
# 5: incentivemonetary - incentivesocial        0.353 -0.076 0.782  2.115   0.069 .  
# 6:     incentiveself - incentivesocial       -0.235 -0.720 0.250 -1.246   0.255    

fit  <- lrm(accuracy_score ~ lateral + incentive, Qualtrics, subset = attentionchecks == 2)
comp <- glht(fit, linfct = c("`lateral=TRUE`=0"))
#              rn coefficients   lwr  upr tstat p.value    
# 1: lateral=TRUE        0.431 0.252 0.61  4.73       0 ***

fit  <- lrm(accuracy_score ~ restraint + incentive, Qualtrics, subset = attentionchecks == 2)
comp <- glht(fit, linfct = c("`restraint=TRUE`=0"))
#                rn coefficients   lwr   upr tstat p.value    
# 1: restraint=TRUE        0.546 0.323 0.768   4.8       0 ***

comparison.summary(comp, "fdr")


# monetary versus non-monetary
fit <- lrm(accuracy_score~incentive, data = Qualtrics[,.(accuracy_score,incentive = factor(incentive, 
                                                   levels = c("control","monetary","self"       ,"social"), 
                                                   labels = c("control","monetary","nonmonetary","nonmonetary")))])
comp <- glht(fit, linfct = c("`incentive=monetary` = 0","`incentive=nonmonetary` = 0","`incentive=monetary` - `incentive=nonmonetary` = 0"))
comparison.summary(comp, "fdr")
#                   coefficients    lwr   upr tstat p.value    
# 1:             monetary  0.495  0.308 0.681  6.21   0.000 ***
# 2:          nonmonetary  0.157 -0.002 0.315  2.31   0.021 *  
# 3: monetary-nonmonetary  0.338  0.176 0.501  4.87   0.000 ***


# mixed effects
fit  <- clmm(ordered(accuracy_score) ~ incentive + (1 | tag), Qualtrics, link = "logit")
#           contrast estimate asymp.LCL asymp.UCL z.ratio p.value    
# 1:   monetary-self    0.315     0.099     0.531    3.85   0.000 ***
# 2: monetary-social    0.399     0.186     0.612    4.94   0.000 ***
# 3:        monetary    0.514     0.301     0.727    6.37   0.000 ***
# 4:     self-social    0.084    -0.126     0.294    1.05   0.291    
# 5:            self    0.199    -0.011     0.408    2.50   0.019 *  
# 6:          social    0.115    -0.092     0.321    1.47   0.171   

fit <- glmer(lateral ~ incentive + (1|tag), Qualtrics, family = "binomial")
#           contrast estimate asymp.LCL asymp.UCL z.ratio p.value    
# 1:   monetary-self    0.559     0.190     0.929   3.992   0.000 ***
# 2: monetary-social    0.486     0.122     0.849   3.526   0.001 ***
# 3:        monetary    0.619     0.243     0.996   4.340   0.000 ***
# 4:     self-social   -0.073    -0.474     0.327  -0.484   0.701    
# 5:            self    0.060    -0.352     0.473   0.385   0.701    
# 6:          social    0.134    -0.274     0.541   0.866   0.580

fit <- glmer(restraint ~ incentive + (1|tag), Qualtrics, family = "binomial")
#           contrast estimate asymp.LCL asymp.UCL z.ratio p.value    
# 1:   monetary-self    0.566     0.101     1.031   3.210   0.004 ** 
# 2: monetary-social    0.351    -0.089     0.792   2.103   0.071 .  
# 3:        monetary    0.612     0.139     1.086   3.416   0.004 ** 
# 4:     self-social   -0.215    -0.710     0.280  -1.145   0.302    
# 5:            self    0.046    -0.478     0.571   0.233   0.816    
# 6:          social    0.261    -0.241     0.764   1.372   0.255

fit <- glmer(correct ~ incentive + (1|tag), Qualtrics, family = "binomial")
# contrast        estimate   asymp.LCL asymp.UCL z.ratio p.value
# monetary           0.410      0.1459     0.674 4.100   <.0001  ***
# self               0.159     -0.1006     0.418 1.610   0.1600 
# social             0.088     -0.1702     0.345 0.900   0.4440 
# monetary-self      0.251     -0.0151     0.518 2.490   0.0260  *
# monetary-social    0.322      0.0572     0.587 3.210   0.0040  **
# self-social        0.071     -0.1890     0.331 0.720   0.4710 

comp <- contrast(emmeans(fit, c("incentive")), list(monetary=c(-1,1,0,0),self=c(-1,0,1,0),social=c(-1,0,0,1),"monetary-self"=c(0,1,-1,0),"monetary-social"=c(0,1,0,-1),"self-social"=c(0,0,1,-1)))
summary(comp, infer = c(T,T), adjust = "fdr")

